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Direct numerical simulations of 
turbulent non-premixed methane-air 
flames modeled with reduced kinetics 

By J. M. Card 1 , J. H* Chen 1 , M. Day 2 , AND S. Mahalingam 3 


Turbulent non-premixed stoichiometric methane-air flames modeled with reduced 
kinetics have been studied using the direct numerical simulation approach. The 
simulations include realistic chemical kinetics, and the molecular transport is mod- 
eled with constant Lewis numbers for individual species. The effect of turbulence 
on the internal flame structure and extinction characteristics of methane-air flames 
is evaluated. Consistent with earlier DNS with simple one-step chemistry, the flame 
is wrinkled and in some regions extinguished by the turbulence, while the turbu- 
lence is weakened in the vicinity of the flame due to a combination of dilatation and 
an increase in kinematic viscosity. Unlike previous results, reignition is observed 
in the present simulations. Lewis number effects are important in determining the 
local stoichiometry of the flame. The results presented in this work are preliminary 
but demonstrate the feasibility of incorporating reduced kinetics for the oxidation 
of methane with direct numerical simulations of homogeneous turbulence to evalu- 
ate the limitations of various levels of reduction in the kinetics and to address the 
formation of thermal and prompt N0 X . 


1* Introduction 

During the 1992 CTR summer program we studied the influence of finite-rate 
chemistry and transient effects on the structure of a turbulent non-premixed flame 
using DNS (Chen et al. 1992a, 1992b, Mahalingam et al. 1994). In this DNS, 
the combustion chemistry was approximated by an Arrhenius single- and two-step 
model. While these simple mechanisms provide useful insights for modeling under 
conditions ranging from near equilibrium to near extinction, they do not, how- 
ever, relate the flame structure and extinction characteristics to the underlying 
elementary chemical-kinetic steps (typically more than 100 reactions) necessary to 
describe the combustion process in hydrocarbons. But the incorporation of such de- 
tailed chemistry for simulations at turbulent Reynolds numbers of practical interest 
is technologically infeasible at the present time. Thus, it is desirable to introduce 
reduced chemical-kinetic mechanisms for the description of the combustion chem- 
istry (Smooke 1991). These mechanisms are deduced from a detailed mechanism by 
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the systematic application of steady-state approximations for intermediate species. 
The combustion kinetics are then described by a few (two to five) global reactions, 
in which only the most elementary kinetic steps are retained and the concentrations 
of various intermediate species are represented by algebraic relations. 

Computations employing reduced mechanisms have successfully been accomplished 
for laminar, planar hydrocarbon, and hydrogen flames (Peters and Rogg 1993). Re- 
cent progress has also been made using reduced kinetics in simulations of reacting 
flows of H 2 - O 2 systems (Baum et al. 1992, and Montgomery et al. 1993). 

The primary objective of this study is to evaluate the effect of different levels of 
reduction of fuel oxidation kinetics on methane-air turbulent diffusion flame struc- 
ture and extinction characteristics. 

2. Reduced chemistry for flame structure 

The oxidation of methane can be described by the following reduced mechanism 
(Seshadri and Peters 1988): 

CH 4 + 2H + H 2 O^CO + 4H 2 , I' 

C0 + H 2 0^C0 2 +H 2 , II' 

2H + M ^ H 2 + M, III' 

0 2 + 3H 2 ^ 2H + 2H 2 0. IV' 

These global reactions were deduced from a 40-some-step starting mechanism by 
the systematic application of steady-state assumptions for the intermediate species. 
The global rates for the above reaction can be expressed in terms of the elementary 


rates shown in Table 1, namely 

WI , = Wl = jfc 38/ [CH 4 ][H], (1) 

w „, = a,„ = (fc 18 //A' 3 )[H]([C0][H 2 0]/[H 2 ] - [C0 2 }/(K 1S /K 3 )), (2) 

win* = wm = MH][0 2 ][M], (3) 

wiv = Jk,/[H]([0 2 ] - [H] 2 [H 2 0] 2 /[H 2 ] 3 A' 1 A' 2 ^D, (4) 


where partial equilibria were assumed for elementary steps (2) and (3). In the above 
rates, the K's axe equilibrium constants, the fc’ s are the elementary rate constants, 
in which the numbering is selected to agree with Peters and Rogg (1993), and where 
t is the concentration of species t. Moreover, the concentration of the third body is 
[ M ] = 6.5[CH 4 ] + 6.5[H 2 0] + 1.5[C0 2 ] + 0.75[CO] -I- 0.4[O 2 ] + 0.4[N 2 ] + 1.0[other]. 
If we further assume that the H atom is in steady state, we get 

0,1V' = WI< + 0,111', 

and the four-step mechanism reduces to the three-step approximation 


CH 4 + 0 2 ^ CO + H 2 + H 2 0, 
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CO + H2O ^ CO2 + H2, II 

0 2 + 2H 2 ^ 2H 2 0, III 

with rates for steps I, II, and II given by Eqs. (1-3), respectively. 

Equations (1),(3), and (4) may be employed to show that 

[H] = (1 - * M/ (CH 4 ]/(*i,[O a l) - k sf [M]/k lf )'' 2 [0al*^IH 1 l*/»(J!r 1 * > ,)»/*jr 1 /(H 8 0]. 

Assuming partial equilibrium for the water-gas shift (eq. II) results in a two- 
step mechanism, and a one-step mechanism may be obtained by assuming H 2 (and 
hence, CO also) to be in steady state. The rate for the global step CH4 + 20 2 
CO2 + 2H2O is taken from Bui-Pham (1992), namely, k = 5.2 x 10 13 exp(— 14906/T) 
(units are moles, cm, s , degK, and KJ/mole). For comparative purposes, future 
calculations will consider the rate constant determined by Puri (1987). 

The four-, three-, and 1-step kinetic mechanisms and the corresponding global 
rates are the basis for the present study. 


3. Governing equations and numerical method 

The governing equations are the continuity, momentum, species, and energy equa- 
tions written in Cartesian tensor notation as follows: 


dp dpuj _ 

Qi *> Q V , 

at oxj 

(5) 

dpui dpuiUj dp drij 

dt dxj dxi dxj ’ 

(6) 

dt + dxj ~ dx’^ pYaVai) + Ua 

(7) 

dpe t d(pe t +p)uj _ d(uj7-y) dq j 

dt dxj dxj dxj 

(8) 

where 

_ / dui duj 2 duk\ 

^\dxj + dxi 3 dxk' ’ 

(9) 

is the stress tensor, 


e < = e + |^“t> e = Y ° ha ~ ~ 

* k = 1 or=X ? 

( 10 ) 


is the total energy per unit mass of the mixture, e is the internal energy per unit 
mass, h Q is the enthalpy of species a given by, 


f T 

h Q =h° a + C pQ {T')dT', 

Jt° 




( 11 ) 
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where h ® and C pa are the corresponding enthalpy of formation and specific heat at 
constant pressure respectively. The heat flux vector is given by 

tYT ^ 

qj = -\ZL + p'£h a Y Q V Qj ( 12 ) 

J a=l 

and the species diffusion velocity is 

V ai = -D aN ±-^, a — 1 ,’”N — 1 (13) 

Y a dxj 

where D a ,v is the binary diffusion coefficient between species a and the iV-th species 
taken to be nitrogen. To ensure that the net diffusion velocity is zero, the constraint 
£a=i = 0 is enforced to obtain V^j- The equation of state is 

p = pRT (14) 

where R is the mixture gas constant given by 

« =*"/«'. < 15 > 


where R° is the universal gas constant, W is the average molecular weight of the 
mixture, and W Q is the species molecular weight. The mixture averaged thermal 
conductivity is modeled through the approximation suggested by Smooke and Gio- 
vangigli (1991), 

A = <?,A(^y C p (T) = jrY a C pa (T ) (16) 


where A = 2.58 x 10 -4 g/cm-sec, r = 0.7, and C p is the specific heat of the gaseous 
mixture. The individual species specific heats are obtained as polynomial functions 
of temperature using the Chemkin thermodynamic database (Kee et al. 1987). The 
diffusion coefficients are obtained by prescription of Lewis numbers for individual 
species through 

= < 17) 

CppLt a 


The Lewis number data is obtained from Smooke and Giovangigli (1991). For the 
one-step mechanism, the Lewis numbers for CH4, O2, CO2, and H2O are 0.97, 1.11, 
1.39, and 0.83, respectively, for the three-step mechanism, the additional Lewis 
numbers for H 2 and CO (0.3 and 1.1, respectively) were employed, and for the 
four-step mechanism, the Lewis number for H is 0.18. The other symbols in these 
equations have the usual meaning. 

For the three-step reduced mechanism, conservation equations for six reacting 
species is considered (CH 4 , O 2 , CO 2 , CO, H 2 , and H 2 O) and the mass fraction of 
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N 2 is obtained through the relationship Y Q = 1, where N = 7. Based on the 

discussion in the previous section, the reaction rates for the six species is computed 
through the following expressions: 

uch 4 = -Wch 4 vi, uq 7 = -Wo 7 [v/ + u >///], lj C o 7 = Wco 7 viu 


and, 


u>co = Wcofoi - a ;//], 


vh 7 = +u;// - 2a;///], 


u;// 2 o = ^ 2 o[w/ - a;// + 2 a;///]. 


(18) 


For the four-step mechanism, N = 8, and seven reacting species (C/f 4 , O2, C0 2 , 
CO, i?2> H 2 0, and H ) are considered. The reaction rate for these species is given 
by: 

UCH 4 = d;o 3 = -W^Oj^/v, 0;c0 2 = W^COa^//, 


o?co = Wco&i — a ;//], a;// 2 - W//J40;/ +w// + a;/// - 3a;/v], (19) 

and, 

^h 2 o = Wh 7 o[2wiv - wii - «/], uh = 2a;/v - 2a;/ — 2a;///. 

The dimensionless form of the governing equations is obtained through definition 
of appropriate reference quantities based on the air stream properties. The sound 
speed at infinity Uoo, corresponding density />oo> and the dynamic pressure pooO^o 
are used for the velocity, density, and pressure reference quantities. The reference 
temperature is (700 — l)Too where 7^ is the ratio of specific heats at temperature 
Too of the air stream. The reference fluid properties used are /ioo, Aoo, and C poo cor- 
responding to air at temperature Too. A reference length scale L te ( is chosen to be 
the distance between the fuel and oxidizer jets in the opposed diffusion flame used to 
initialize all the dependent variables in the computational domain. The equations 
are solved using a standard sixth-order accurate compact finite differencing scheme 
(Lele, 1992) for approximating spatial derivatives, and a third-order Runge-Kutta 
scheme for time advancement. A modified version of the Navier-Stokes Character- 
istic Boundary Condition (NSCBC) procedure originally developed by Poinsot and 
Lele (1992), and suitably modified to account for variable specific heats is imple- 
mented. The boundary conditions are periodic in the y direction and non-reflecting 
in the x direction (see Chen ef a/., 1992a). 

The turbulence field is prescribed by an initial two-dimensional turbulent kinetic 
energy spectrum function 

< 2 °> 

where k is the wavenumber, k$ is the wavenumber corresponding to the most ener- 
getic eddies, and uq is the rms velocity. 
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4. Initial conditions and flame parameters 

The reacting flow field in a steady opposed jet diffusion flame configuration in- 
cluding a full 40 step chemical kinetic mechanism was obtained using the Sandia 
code OPPDIF. Pure methane and air at 1 atmosphere, 300°K with a strain rate 
of approximately 20 sec -1 are the conditions imposed at the fuel and oxidizer jet 
separated by a distance of 1 cm. These correspond to low strain rate conditions. 
The strained laminar flame solution from OPPDIF was then allowed to relax to an 
unstrained condition in the DNS code. The unstrained laminar flame profile was 
used to initialize all of the species mass fractions, temperature, and density fields 
in the turbulence simulations. 

The DNS was initialized in a two-dimensional computational domain with a plane 
laminar diffusion flame in the center of the domain. The fuel stream is on the left 
half of the domain while the oxidizer stream is on the right half of the domain. The 
velocity field was initialized with the turbulent kinetic energy spectrum given in 
(20). The turbulence Reynolds number was taken to be 57 based on the turbulence 
microscale. The flame is thicker than the turbulence microscale but smaller than 
the turbulence integral scale. Two-dimensional simulations were performed on a 
129 by 129 grid. 

5. Results and discussion 

5.1 Laminar flamelei results 

To initialize the direct simulations, the structure for laminar counterflow flamelets 
were calculated for low strain rates (< 100 sec *). The details for such calculations 
may be found in Smooke (1991). Shown in Figs. 1 and 2 are the flamelet structure 
in mixture fraction space for a global one-step reaction and a three-step mechanism, 
respectively. 

A distinctive qualitative difference between the one-step and three-step mech- 
anism is the broadness of the reaction zone. In Fig. lc, it can be seen that the 
consumption of the fuel extends over a much wider range of mixture fraction for 
the one-step mechanism. The thinness of the fuel-consumption layer in multi-step 
mechanisms is one of the primary causes of resolution difficulties in direct sim- 
ulations of turbulent combustion flames. In addition, it can readily be seen in 
Figs, lb and 2b that the fuel and oxygen profiles are quite different. In the one-step 
approximation, fuel leaks through the reaction zone. However, for the three-step 
mechanism, only oxygen leaks through the reaction layer, which is in agreement 
with experiments and numerical calculations employing full detailed mechanisms. 

Preliminary calculations using the four-step mechanism indicate that the qual- 
itative features of the flamelet agree with previous studies. The major expected 
difference between the three- and four-step mechanisms is the H atom profile. In 
the three-step mechanism, all the H is produced on the fuel-lean side of the flamelet 
and is consumed entirely in the fuel-consumption layer. For the four-step approxi- 
mation, the H atom profile is broader, in which the fuel-consumption layer is within 
the H nonequilibrium layer. The monograph by Smooke (1991) may be consulted for 
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Mixture Fraction 

FIGURE 1 . One-dimensional unstrained laminar flame for methane-air one-step 
mechanism: a) temperature, b) species mass fractions, c) species reaction rates. 
( H 2 0, C02, — 02 ? — CH 4 ). 

more detailed discussion on the structure of laminar methane flames using reduced 
chemistries. 

5.2 Turbulent flame structure 

Preliminary results are presented here for the DNS using the global one-step mech- 
anism described in section 2. Simulations are currently underway for the three- and 
four- step mechanisms. In order to resolve the fuel consumption layer the resolution 
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FIGURE 2. One- dimensional strained laminar flame for methane-air three-step 

mechanism: a) temperature, b) species mass fractions, c) species reaction rates. 

( H 2 0, H 2 * 40, CO, C0 2 , 0 2 , CH a ). 

requirements are much greater for the three- and four-step mechanisms than for the 
one-step mechanism. 

The turbulent two-dimensional DNS was postprocessed at 1.0 and 1.5 turbulent 
eddy turn times. This is sufficient time for the initial adjustment to occur between 
the imposed turbulence spectrum and the laminar flame. Instantaneous images of 
the vorticity magnitude, temperature, reaction rate, and species concentrations at 
t = 1.0 eddy turn time are shown in Fig. 3. It is evident from Fig. 3a that the 
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H I ! i “ 1 I I ! x 

FIGURE 4. CH 4 reaction rate at 1.5 eddy turn times. 

vorticity is greatly dampened in the vicinity of the reaction zone. The turbulence 
was initialized at all locations in the computational domain except near the outflow 
boundaries normal to the flame. Due to a combination of dilatational effects and 
a substantial increase in the kinematic viscosity with temperature, vorticity in the 
vicinity of the flame is dampened. Note in Fig. 3c that the vorticity is dampened 
more on the fuel side than on the oxidizer side. This is due to the lower density 
(higher temperature) that exists on the rich side of the flame which results in a 
higher kinematic viscosity on that side. 

While the effect of the heat release is to weaken the turbulence, the effect of 
the turbulence is to wrinkle the flame. There Me regions where the flame is locally 
extinguished as shown by the discontinuities in the reaction rate contour lines shown 
in Fig. 3c. At a later time in the simulation, t = 1.5 eddy turn times, the conditions 
are such that locally the strain rate is reduced in regions with hot products due to 
the heat release. The hot products act as an ignition source and in locations where 
fuel and oxidizer are present, the flame starts to burn once again. A comparison of 
Fig. 4 and Fig. 3c shows a region where reignition occurs. Previous simulations with 
simple chemistry (Chen 1992a, 1992b) did not report the occurrence of reignition. 

It is also observed from the reaction rate contours and the stoichiometric mixture 
fraction line (Fig. 3c) that the flame does not always burn at the stoichiometric 
mixture fraction. Near the peak reaction rates the flame is stoichiometric; however, 
in regions that are curved toward the oxidizer stream the flame tends to burn lean, 
whereas in regions that are curved toward the fuel stream the flame tends to burn 
rich. We suspect that this is a Lewis number effect and that it is particularly pro- 
nounced in regions having the largest curvature. In regions that are curved toward 
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Figure 5. Scatter plots at 1.0 eddy turn time for one-step mechanism: 
a) temperature, b) CH 4 reaction rate, c) CH 4 and O 2 mass fraction (□ x 0 2 ), 

d) C02 mass fraction, and e) H^O mass fraction. 
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the oxidizer, the defocusing of methane overcomes the greater mass diffusivity of 
methane ( Le = 0.97) over oxygen ( Le = 1.11). Hence, the flame tends to burn 
lean. On the other hand, in regions that are curved towards the fuel, the additive 
effect of focusing the methane combined with its larger mass diffusivity leads to a 
situation where there is an abundance of methane at the reaction zone. Hence, in 
this situation the flame tends to burn rich, and eventually, as the curvature becomes 
too large, the flame becomes too rich to burn, at which point local extinction occurs 
first at the flame tip. Since the reactant Lewis numbers in the DNS can readily be 
interchanged, these influences will be studied in greater detail in future simulations. 

The reactant and product concentrations are presented in Figs. 3d-3f. Note that 
CH 4 is more diffusive them 02- The products CO 2 and H 2 O follow the temperature 
(Fig. 3b) quite closely; hence, for the sake of brevity, only CO 2 is presented in Fig. 3f. 

Scatter plots of the temperature, CH 4 reaction rate, and species concentrations 
in mixture fraction space are shown in Fig. 5. The mixture fraction, f, is defined 
to be zero in the oxidizer stream and unity in the fuel stream. The maximum 
temperature occurs near the stoichiometric mixture fraction, / = 0.057. The CH 4 
reaction rate peak is shifted slightly toward the fuel side. Consistent with the 
physical representation of the flame, the presence of both high and low values of the 
temperature and reaction rate in Figs. 5a and 5b suggests that there are regions in 
the flame that are fully burning as well as other regions that are being extinguished 
by the local strain-rate field. From Fig. 5c, note that there is leakage of both 
CH 4 and 0 2 in this one-step model. Also, a significant portion of the flame is 
undergoing various stages of extinction as indicated by the spread of points between 
the equilibrium and frozen flow limits. 

6. Concluding remarks 

Preliminary results from DNS demonstrate that two-dimensional simulations of 
turbulent non-premixed methane- air flames modeled with reduced chemistry are 
feasible. The global one-step simulation required 50 cpu hours on a Cray-YMP. Ex- 
tensions to three- and four- step mechanisms requires resolving the fuel consumption 
layer of the flame, a region that is approximately an order of magnitude smaller 
than the overall flame thickness. Simulations on a uniform grid will have enormous 
resolution requirements if the internal flame structure is to be fully resolved. This 
suggests a real need for a local adaptive mesh refinement capability for simulations 
of hydrocarbon flames. 

Once the three- and four-step simulations are completed, detailed comparisons 
of the internal flame structure and extinction characteristics will be made. In par- 
ticular, strain-rate and curvature statistics will be obtained and comparisons of 
the turbulence simulations with steady one-dimensional strained laminar flame pre- 
dictions will be made. Finally, reduced mechanisms for the formation of thermal 
and prompt NO x will be evaluated from the major species concentrations and the 
temperature. 
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